library(car)
library(lme4)
library(lmerTest)
library(ggplot2)
library(dplyr)
library(lattice)
set.seed(1234)
Blackmore$log.exercise <- log(Blackmore$exercise + 5/60, 2)
Blackmore$age0 <- Blackmore$age - 8
bm.lme.1 <- lmer(log.exercise ~ age0*group + (age0|subject),
                 data=Blackmore)
summary(bm.lme.1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log.exercise ~ age0 * group + (age0 | subject)
##    Data: Blackmore
## 
## REML criterion at convergence: 3614.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7349 -0.4245  0.1228  0.5280  2.6362 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  subject  (Intercept) 2.08384  1.4436        
##           age0        0.02716  0.1648   -0.28
##  Residual             1.54775  1.2441        
## Number of obs: 945, groups:  subject, 231
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        -0.27602    0.18237 236.27186  -1.514   0.1315    
## age0                0.06402    0.03136 237.54226   2.041   0.0423 *  
## grouppatient       -0.35400    0.23529 233.73375  -1.505   0.1338    
## age0:grouppatient   0.23986    0.03941 221.23037   6.087 5.03e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) age0   grpptn
## age0        -0.489              
## grouppatint -0.775  0.379       
## ag0:grpptnt  0.389 -0.796 -0.489

Leverage

Leverage is a measure of how much an individual influenced the model fit. It is a one-number summary of how different the model fit would be if the given individual was excluded, compared to the model fit where the individual is included.

The hat-value, \(h_i\), is a common measure of leverage and they can be obtained using the function hatvalues(modelObjectName). The response-variable values are not at all involved in determining leverage.

summary(hatvalues(bm.lme.1))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1724  0.2170  0.2672  0.2834  0.3461  0.4629

Outliers

Discrepant observations usually have large residuals. The outlierTest() function in the car package performs a Bonferroni adjustment to the p value for testing the statistical significance of the studentized residuals. The studentized residual follows a \(t\)-distribution with \(n-k-2\) degrees of freedom. Studentized residuals can be obtained using the function rstudent(modelObjectName).

summary(rstudent(bm.lme.1))
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -3.466640 -0.516408  0.144848  0.005752  0.631648  3.263598
outlierTest(bm.lme.1)  # function from the car package
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 848 -3.46664         0.00055097      0.52067

A graphical approach to detecting outliers is to construct a QQ plot for the studentized residuals. QQ plots are useful for examining the distribution of the residuals, but are also effective in displaying the tail behavior of the residuals: outliers, skewness, heavy tails, or light tails all show up clearly.

qqmath(bm.lme.1)             # function from the lattice package

Influence

Influence on the regression coefficients combines leverage and discrepancy. The following heuristic formula helps to distinguish among the three concepts of influence, leverage, and discrepancy (‘outlyingness’):
Influence on Coefficients = Leverage × Discrepancy

Plot the studentized residuals against the hat values and look for observations for which both are large.

rstudent <- rstudent(bm.lme.1)
leverage <- hatvalues(bm.lme.1)
ggplot( , aes(leverage, rstudent)) +
  geom_vline(size = 2, colour = "white", xintercept = 0) +
  geom_hline(size = 2, colour = "white", yintercept = 0) +
  geom_point() + 
  ggtitle("Studentized Residual vs. Leverage")

To obtain influence measures for models fit using the lme4 package, use the influence function and save them in an object. Note that this takes a little time to run so you may want to use cache=TRUE.

infMeasures <- influence(bm.lme.1, "subject")

You can then use the influenceIndexPlot() function from the car package to generate diagnostic plots.

influenceIndexPlot(infMeasures)

The bottom plot is for Cook’s Distance. The top 4 plots are for the DFBETAS.

The DFBETAS is most direct measure of influence and simply expresses the impact on each coefficient of deleting each individual in turn. Large values of DFBETAS indicate individuals that are influential in estimating a given parameter. One problem associated with using DFBETAS is there are a large number of them (one for each individual and each parameter estimated). It is useful to have a single summary index of influence.

Cook’s Distance is a single summary index of influence produced for each individual and can be obtained using the function cooks.distance().

Finally, specifying vars="var.cov.comps" in the influenceIndexPlot() function will give diagnostic plots for the variance covariance components. In our example, the random intercept variance, random slope variance, the covariance between the random intercepts and slopes, and the error variance.

influenceIndexPlot(infMeasures, vars="var.cov.comps")

An ideal case

The plot below is a great example of a Studentized Residuals vs Leverage plot in which we see no evidence of outliers. None of the points come close to having both high residual and leverage.

## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

influenceIndexPlot(lm.good)

Assumptions of the Mixed-Effects Model

\(\epsilon_{ij} \sim N(0,\sigma_{\epsilon}^2)\)

\(\sigma^2_{b_0}\) and \(\sigma^2_{b_1}\) are allowed to covary but both are independent of \(\epsilon_{ij}\). Thus, \(\mathbf{b}_i \sim \mathbf{N}_q(\mathbf{0,G})\).

Are the level-one residuals normally distributed?

The Normal Q-Q plot is a probability plot of the studentized residuals against the values that would be expected under normality. If the normality assumption is met, the points on this graph should fall on a straight 45-degree line.

qqPlot(rstudent(bm.lme.1))  # function from the car package

## 848 922 
## 489 563

An ideal case

Below is QQ plot in which the residuals satisify the normality assumption.

## [1] 111 414

The residuals here are a perfect match to the diagonal line, and thus are normally distributed. This is an example of a Normal QQ plot that is as perfect as it gets.

Lighter Tails

In the next example, we see a QQ plot where the residuals deviate from the diagonal line in both the upper and lower tail. This plot indicates that the tails are ‘lighter’ (have smaller values) than what we would expect under the normal distribution. This is indicated by the points forming a “flatter” line than than the diagonal.

## [1]  53 695

Heavier Tails

In this final example, we see a QQ plot where the residuals again deviate from the diagonal line in both the upper and lower tail. Unlike the previous plot, in this case we see that the tails are observed to be ‘heavier’ (have larger values) than what we would expect under the normal distribution. This is indicated by the points forming a “steeper” line than the diagonal.

## [1] 795  37

Do the level-one residuals have constant variance?

Non-constant error variance is sometimes termed ‘heteroscedasticity.’ Non-constant error variance is a serious problem only when it is relatively extreme. Furthermore, unlike the ordinary regression model, the mixed-effects model can accommodate heteroscedasticity and dependencies among the errors if we do not assume that \(\mathbf{R}_i = \sigma^2_{\epsilon} \mathbf{I}_{n_i}\)

It is common for error variance to increase as the expectation of \(Y\) grows larger, or there may be a systematic relationship between error variance and a particular \(X\). The former situation can often be detected by plotting studentized residuals against fitted values; the latter by plotting residuals against each \(X\).

plot(bm.lme.1)

Scale - Location

It often helps to plot \(\sqrt{|\widehat{\epsilon}_i|}\) against \(\widehat{Y}\), which is referred to as a scale-location plot.

plot(bm.lme.1, sqrt(abs(resid(.))) ~ fitted(.))

Example in which the constant error variance assumption does not hold

Here is a scatterplot of the data.

Here is what the Studentized Residual vs. Fitted plot looks like in this case.

The distribution of the residuals is quite well concentrated around 0 for small fitted values, but they get more and more spread out as the fitted values increase. This is an instance of “increasing variance”.

Here is the scale-location plot, which tells a similar story.

Are the level-two random coefficients normally distributed with mean 0?

randEffects <- ranef(bm.lme.1)
plot(randEffects)
## $subject

Using the QQ plot:

qqmath(randEffects)
## $subject

Linearity

Next, we must assess whether the predictor variables are linearly related to the outcome. If the outcome variable is linearly related to the predictor variables, there should be no systematic relationship between the residuals and the model predicted (that is, fitted) values.

The linearity assumption is easily checked by plotting the studentized residuals vs. the fitted values. If linearity holds, there is no systematic relationship. But if linearity does not hold, it will be apparent.

plot(bm.lme.1)

Here is a scatterplot of two variables, \(X\) and \(Y\), which are not linearly related.

Suppose we fit a regression of \(Y\) on \(X\) and obtain the fitted values and studentized residuals. Here is a plot of the studentized residuals vs. the fitted values. It is clear that there is a systematic relationship.

Plotting studentized residuals against each \(X\) may also be helpful for detecting departures from linearity.

The plot suggests a polynomial relationship so let’s include a quadratic term in the regression.

lm.poly <- lm(y.poly ~ poly(x, 2))

rstudent <- rstudent(lm.poly)
fitted <- fitted(lm.poly)

ggplot( , aes(fitted, rstudent)) +
  geom_point() +
  geom_smooth(se=FALSE) +
  xlab("Fitted values") + ylab("Studentized Residuals") +
  ggtitle("Studentized Residual vs Fitted Plot")

There is no longer a systematic relationship.

Nonlinearities can often be corrected by considering polynomial or spline functions.

Summary

Non-normality can often be detected by examining the distribution of the studentized residuals, and frequently can be corrected by transforming the data.

The following four plots are important diagnostic tools in assessing whether the model assumptions hold and in detecting influential observations.

Studentized Residuals vs. Fitted Values

When a linear mixed-effects model is appropriate, we expect

  1. the studentized residuals will have constant variance when plotted against fitted values; and

  2. the studentized residuals and fitted values will be uncorrelated.

If there are clear trends in this plot, or the plot looks like a funnel, these are clear indicators that the given model is inappropriate.

Normal QQ plot

The studentized residuals should be approximately normally distributed. The QQ plot is also useful for detecting outliers.

Scale-location plot

This is another version of the studentized residuals vs fitted values plot. There should be no discernible trends in this plot.

Studentized Residuals vs Leverage

Points with a high residual (poorly described by the model) and high leverage (high influence on model fit) are potentially influential.